******************
*Code for replication of tables and figures in "Exercise Imporves Academic Performance"
*Note1: Variables with _b suffix are from the baseline survey; variables with _m suffix are from the midterm survey.
*Note2: This code is intended as a step-by-step procedure, and not intended to run from start to beginning recreating all files.
*Note3: You may need to install additional packages to perfrom some of the analysis
*Note3: For details on the variables please consult the paper and the full surveys published in the online appendix.
******************

*******************************************************************************************************************************************
*IMPORTANT: To comply with requirements to preserve the anonymity of participants we had to remove the variable "age" (which was used as a control variable in the paper) from the publicly available data set. The regression results with control variables on the publicly available dataset can therefore vary slightly from the results reported in the paper. As shown in Table 1 in the readme file, however, removing the age variable from the regressions with control variables has no consequences for our main treatment effects reported in Table 3, 4, 6 and 7 of the paper. The point estimates are almost identical, and the significance levels remain unchanged.     
*********************************************************************************************************************************************

clear all 
set more off
set scheme s1mono, permanent
global pm = char(177)

use gym_JPE_final_anonymized.dta


***std Gym
sum visitspersem visited

gen std_gymvisits=visitspersem/9.044197
gen std_visited=visited/.4916769

gen std_gymvisits_sq=std_gymvisits^2

***std Academic
sum CSP GA
gen std_CSP=CSP/11.73987
gen std_GA=GA/.9673969

sum CSP_retake CSP_pf 
gen std_CSP_retake=CSP_retake/5.933093  
gen std_CSP_pf=CSP_pf/5.903643 

sum RSP ESP
gen std_RSP=RSP/14.81521
gen std_ESP=ESP/11.88121

sum CSP_F15 GA_F15 RSP_F15

gen std_CSP_F15=CSP_F15/10.63542
gen std_GA_F15=GA_F15/.9340003
gen std_RSP_F15=RSP_F15/12.21625

***Generate and standardize intermediate outcomes 
gen sum_alcohol_m=vine_m+beer_m+otheralcohol_m 
egen std_alcohol_m=std(sum_alcohol_m)

egen std_bedaftermidnight_m=std(bedaftermidnight_m) 
egen std_daystired_m=std(daystired_m)
egen std_happy_health_m=std(happy_health_m)

gen health_index_m=std_happy_health_m-std_daystired_m-std_bedaftermidnight_m
gen health_index_m_alc=std_happy_health_m-std_daystired_m-std_bedaftermidnight_m-std_alcohol_m

gen dontsticktoplans_m=11-sticktoplans_m
gen selfcontrol_m=11-(acttoofast_m+difficultresisttemptations_m+procrastinate_m+dontsticktoplans_m)/4

sum health_index_m health_index_m_alc hoursstudied_m happy_general_m selfcontrol_m timesexercised_m

gen std_lifestyle=health_index_m/2.058044 
gen std_lifestyle_alc=health_index_m_alc/2.374391 
gen std_study=hoursstudied_m/17.35442
gen std_happy=happy_general_m/1.944128 
gen std_selfcontrol=selfcontrol_m/1.478326
gen std_exercise=timesexercised_m/2.454513

***Above median on at least one dimension 
gen combined_above=above_median_health+above_median_study+above_median_happiness+above_median_selfcontrol

gen above=1 if inrange(combined_above,1,4)
replace above=0 if combined_above==0

***Above median pre-exercise
gen above_exercise=.
replace above_exercise=0 if inrange(timesexercised_b, 0,1)
replace above_exercise=1 if inrange(timesexercised_b, 2,14)

***Interactions
gen t_above_health=treated*above_median_health
gen t_above_happy=treated*above_median_happiness
gen t_above_study=treated*above_median_study
gen t_above_self=treated*above_median_selfcontrol
gen t_above_exercise=treated*above_exercise

gen t_above=treated*above

***attrition variable
gen allbackground=1 
replace allbackground=0 if above_median_study==. | above_median_health==. | above_median_happiness==. | above_median_selfcontrol==.

gen alloutcome=1 
replace alloutcome=0 if std_lifestyle_alc==. | std_selfcontrol==. | std_happy==. | std_study==.


********************
*Tables in paper and appendix
********************


********
***Paper
********

***balance test (Table 2) 
*(requires orth_out program "ssc install orth_out")
orth_out female_b yearofstudy_b UIB ///
selfcontrol_b health_index hoursstudied_b happy_general_b  RSP ///
using sum_stat_short.tex,by(treated) test overall count latex full replace

***Main Gym (Table 3)
eststo clear
eststo: quietly reg std_gymvisits treated,r
eststo: quietly reg std_gymvisits treated female_b yearofstudy_b UIB,r
eststo: quietly reg std_gymvisits treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg visited treated,r
eststo: quietly reg visited treated female_b yearofstudy_b UIB,r
eststo: quietly reg visited treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
esttab , se star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("Visits" "Visits" "Visits" "Visit (1/0)" "Visit (1/0)" "Visit (1/0)") title(Regressions: Gym (std)) replace

*add p-values from MHT corrections
*RW-method (requires rwolf program "ssc install rwolf")
rwolf std_gymvisits visited, indepvar(treated) controls(female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol) reps(10000) seed(19283746)

*HB (Holm from output) and LSV-method (Thm3_1 from output). NOTE: code provided by Atom Vayalinkal (the mhtexp2 package is available at https://github.com/vayalinkal/mhtexp2, we used a demo version of it).
*(requires moremata program "ssc install moremata"")
preserve

mata: mata mlib index

// drop obs that don't have data for controls or subgroup identifiers
local allvars female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol
foreach v in `allvars'{
                drop if `v' == .

}

disp _N
local sN = _N
disp `sN'
local B = 3000
local stu = 1
mata: rseed(0)
// pregenerate randomization matrix to save time
mata: idbootmat = runiformint(`sN', `B', 1, `sN')

mhtexp2 std_gymvisits visited, treatment(treated) controls(female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol) 


***Main Academic (Table 4)
eststo clear
eststo: quietly reg std_CSP treated,r
eststo: quietly reg std_CSP treated female_b yearofstudy_b UIB std_RSP,r
eststo: quietly reg std_CSP treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP,r
eststo: quietly reg std_GA treated,r
eststo: quietly reg std_GA treated female_b yearofstudy_b UIB std_RSP,r
eststo: quietly reg std_GA treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP,r
esttab , se star(* 0.1 ** 0.05 *** 0.01)  nonumbers mtitles("CSP" "CSP" "CSP" "GA" "GA" "GA") title(Regressions: Academic (std)) replace

*add p-values from MHT corrections
*RW-method 
rwolf std_CSP std_GA, indepvar(treated) controls(female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP) reps(10000) seed(19283746)

*HB (Holm from output) and LSV-method (Thm3_1 from output): CSP Outcome 
preserve

mata: mata mlib index

// drop obs that don't have data for controls or subgroup identifiers
local allvars female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP
foreach v in `allvars'{
                drop if `v' == .

}

disp _N
local sN = _N
disp `sN'
local B = 3000
local stu = 1
mata: rseed(0)
// pregenerate randomization matrix to save time
mata: idbootmat = runiformint(`sN', `B', 1, `sN')

mhtexp2 std_CSP std_GA, treatment(treated) controls(female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP) 

*HB and LSV-method: GA outcome
preserve

keep if std_GA!=.

mata: mata mlib index

// drop obs that don't have data for controls or subgroup identifiers
local allvars female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP
foreach v in `allvars'{
                drop if `v' == .

}

disp _N
local sN = _N
disp `sN'
local B = 3000
local stu = 1
mata: rseed(0)
// pregenerate randomization matrix to save time
mata: idbootmat = runiformint(`sN', `B', 1, `sN')

mhtexp2 std_CSP std_GA, treatment(treated) controls(female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP) 


*** CSP and Gym (Table 5)
eststo clear
eststo: quietly reg std_CSP std_gymvisits std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_CSP std_gymvisits std_gymvisits_sq std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_CSP std_gymvisits std_gymvisits_sq treated std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_CSP std_gymvisits std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol if treated==1,r
eststo: quietly reg std_CSP std_gymvisits std_gymvisits_sq std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol if treated==1,r

esttab , se star(* 0.1 ** 0.05 *** 0.01)  nonumbers mtitles("std CSP" "std CSP" "std CSP" "std CSP: Treated" "std CSP: Treated") drop(std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol)  title(Controls: std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol)

*** Intermediary (Table 6)
eststo clear
eststo: quietly reg std_lifestyle_alc treated,r
eststo: quietly reg std_lifestyle_alc treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_lifestyle_alc treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_gymvisits std_gymvisits_sq,r
eststo: quietly reg std_selfcontrol treated,r
eststo: quietly reg std_selfcontrol treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_selfcontrol treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_gymvisits std_gymvisits_sq,r
eststo: quietly reg std_happy treated,r
eststo: quietly reg std_happy treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_happy treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_gymvisits std_gymvisits_sq,r
eststo: quietly reg std_study treated,r
eststo: quietly reg std_study treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_study treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_gymvisits std_gymvisits_sq,r
esttab , se star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("Lifestyle" "Lifestyle" "Lifestyle""Selfcontrol" "Selfcontrol" "Selfcontrol" "Happy" "Happy" "Happy" "Study" "Study" "Study") title(Regressions other outcomes (std):) replace

*MHT rwolf
rwolf std_lifestyle_alc std_selfcontrol std_happy std_study, indepvar(treated) controls(female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol) reps(10000) seed(19283746)
*MHT Holm and LSV-method
preserve

keep if alloutcome==1

mata: mata mlib index

// drop obs that don't have data for controls or subgroup identifiers
local allvars female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol
foreach v in `allvars'{
                drop if `v' == .

}

disp _N
local sN = _N
disp `sN'
local B = 3000
local stu = 1
mata: rseed(0)
// pregenerate randomization matrix to save time
mata: idbootmat = runiformint(`sN', `B', 1, `sN')

mhtexp2 std_lifestyle_alc std_selfcontrol std_happy std_study, treatment(treated) controls(female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol) 

*bounds (requires leebounds program "ssc install leebounds")
preserve

leebounds std_lifestyle_alc treated, cie vce(bootstrap, reps(3000) nodots) level(95)
leebounds std_selfcontrol treated, cie vce(bootstrap, reps(3000) nodots) level(95)
leebounds std_happy treated, cie vce(bootstrap, reps(3000) nodots) level(95)
leebounds std_study treated, cie vce(bootstrap, reps(3000) nodots) level(95)

***Hetro Academic (Table 7)
eststo clear
eststo: quietly reg std_CSP treated above t_above female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_CSP treated t_above_health above_median_study  above_median_health  above_median_happiness  above_median_selfcontrol female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above_health
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_CSP treated t_above_self above_median_study  above_median_health  above_median_happiness  above_median_selfcontrol female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above_self
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_CSP treated t_above_happy above_median_study  above_median_health  above_median_happiness  above_median_selfcontrol female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above_happy
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_CSP treated t_above_study above_median_study  above_median_health  above_median_happiness  above_median_selfcontrol female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above_study
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_GA treated above t_above female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)

esttab , se stat(skal skalse N) star(* 0.1 ** 0.05 *** 0.01) drop(female_b yearofstudy_b UIB std_RSP) nonumbers mtitles("CSP" "CSP" "CSP" "CSP" "CSP" "GA") title(Regressions: heterogeneity (std)) replace

**************
*Appendix
**************

***Correlation table (Table A1)
gen sum_alcohol_b=vine_b+beer_b+otheralcohol_b
gen thinkbefore_b=11-acttoofast_b
gen resist_b=11-difficultresisttemptations_b
gen dontprocrastinate_b=11-procrastinate_b

eststo clear
estpost corr bedaftermidnight_b daystired_b sum_alcohol_b happy_health_b thinkbefore_b resist_b dontprocrastinate_b sticktoplans_b, matrix casewise
eststo correlation
esttab correlation using correlations.csv, star(* 0.1 ** 0.05 *** 0.01)  replace not unstack compress
esttab .,  star(* 0.1 ** 0.05 *** 0.01)  replace not unstack compress

***long sum stat (Table A2)
orth_out female_b yearofstudy_b UIB ///
selfcontrol_b health_index hoursstudied_b happy_general_b  RSP  ///
work_b hoursworked_b memberofothergym_b gymmemberbefore_b timesexercised_b VPS_F15 CSP_F15 GA_F15 ///
using sum_stat_long.tex,by(treated) test overall count latex full replace

***Gym effects by barrier to exercise (Table A3)
gen  expensivebuttime_8=1 if inrange(SIBcardexpensive_b,8,10)
replace expensivebuttime_8=0 if inrange(SIBcardexpensive_b,1,7)
replace expensivebuttime_8=0 if inrange(notimetoexercise_b, 8,10)

*code NA's as zero (does not impact result but maintain full sample)
gen expensivebuttime_8_na=expensivebuttime_8
replace expensivebuttime_8_na=0 if expensivebuttime_8_na==.

gen treated_expensivebuttime_8=treated*expensivebuttime_8
gen treated_expensivebuttime_8_na=treated*expensivebuttime_8_na

gen  notime_8=1 if inrange(notimetoexercise_b,8,10)
replace notime_8=0 if inrange(notimetoexercise_b,1,7)

*code NA's as zero (does not impact result but maintain full sample)
gen notime_8_na=notime_8
replace notime_8_na=0 if notime_8_na==.

gen treated_notime_8=treated*notime_8
gen treated_notime_8_na=treated*notime_8_na

eststo clear
eststo:  quietly reg std_gymvisits treated expensivebuttime_8_na treated_expensivebuttime_8_na female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol, r
lincom treated+treated_expensivebuttime_8_na
estadd scalar skal_cost r(estimate)
estadd scalar skalse_cost r(se)
eststo:  quietly reg std_gymvisits treated notime_8_na treated_notime_8_na female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol, r
lincom treated+treated_notime_8_na
estadd scalar skal_time r(estimate)
estadd scalar skalse_time r(se)
eststo:  quietly reg std_gymvisits treated expensivebuttime_8_na treated_expensivebuttime_8_na notime_8_na treated_notime_8_na female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol, r
lincom treated+treated_expensivebuttime_8_na
estadd scalar skal_cost r(estimate)
estadd scalar skalse_cost r(se)
lincom treated+treated_notime_8_na
estadd scalar skal_time r(estimate)
estadd scalar skalse_time r(se)
esttab,se stat(skal_cost skalse_cost skal_time skalse_time N)  star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("std GYM" "std GYM" "std GYM w/cont" "std GYM" "std GYM" "std GYM w/cont") drop (female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol) nogaps

*robustness (w/o NA's)
eststo clear
eststo:  quietly reg std_gymvisits treated expensivebuttime_8 treated_expensivebuttime_8 female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol, r
lincom treated+treated_expensivebuttime_8
estadd scalar skal_cost r(estimate)
estadd scalar skalse_cost r(se)
eststo:  quietly reg std_gymvisits treated notime_8 treated_notime_8 female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol, r
lincom treated+treated_notime_8
estadd scalar skal_time r(estimate)
estadd scalar skalse_time r(se)
eststo:  quietly reg std_gymvisits treated expensivebuttime_8 treated_expensivebuttime_8 notime_8 treated_notime_8_na female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol, r
lincom treated+treated_expensivebuttime_8
estadd scalar skal_cost r(estimate)
estadd scalar skalse_cost r(se)
lincom treated+treated_notime_8
estadd scalar skal_time r(estimate)
estadd scalar skalse_time r(se)
esttab,se stat(skal_cost skalse_cost skal_time skalse_time N)  star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("std GYM" "std GYM" "std GYM w/cont" "std GYM" "std GYM" "std GYM w/cont") drop (female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol) nogaps

*** Academic perfromance controlling for Fall 2015 (Table A4)

eststo clear
eststo: quietly reg std_CSP treated std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_CSP treated std_CSP_F15 std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_CSP treated std_CSP_F15 std_GA_F15 std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_GA treated std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_GA treated std_GA_F15 std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_GA treated std_GA_F15 std_CSP_F15 std_RSP female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r

esttab , se star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles( "CSP" "CSP" "CSP" "GA"  "GA" "GA") title(Regressions: Control for Fall 2015 (std)) drop(female_b yearofstudy_b UIB above_median_health above_median_study above_median_happiness above_median_selfcontrol std_RSP) replace

*** Academic perfromance including alchol in lifestyle index (Table A5)
eststo clear
eststo: quietly reg std_CSP treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP,r
eststo: quietly reg std_CSP treated female_b yearofstudy_b UIB above_median_study above_median_health2 above_median_happiness above_median_selfcontrol std_RSP,r
eststo: quietly reg std_GA treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP,r
eststo: quietly reg std_GA treated female_b yearofstudy_b UIB above_median_study above_median_health2 above_median_happiness above_median_selfcontrol std_RSP,r
esttab , se star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("CSP w/o" "CSP w/" "GA w/o" "GA w/") drop(female_b yearofstudy_b UIB above_median_health above_median_health2 above_median_study above_median_happiness above_median_selfcontrol std_RSP) title(Regressions: incl alcohol in health) replace


*** Subparts lifestyle (Table A6)
sum bedaftermidnight_m daystired_m happy_health_m sum_alcohol_m

gen std_bed=bedaftermidnight_m/2.152701
gen std_tired=daystired_m/1.869221
gen std_health=happy_health_m/2.172137 
gen std_alcohol=sum_alcohol_m/6.490323 

*above median subindex
tabstat bedaftermidnight_b daystired_b sum_alcohol_b happy_health_b,stat(p50)

gen am_bed=0
replace am_bed=1 if bedaftermidnight_b>3
replace am_bed=. if bedaftermidnight_b==.

gen am_tired=0
replace am_tired=1 if daystired_b>3
replace am_tired=. if daystired_b==.

gen am_alc=0
replace am_alc=1 if sum_alcohol_b>4
replace am_alc=. if sum_alcohol_b==.

gen am_health=0
replace am_health=1 if happy_health_b>7
replace am_health=. if happy_health_b==.

eststo clear
eststo: quietly reg std_lifestyle_alc treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_bed treated am_bed female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_tired treated am_tired female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_alcohol treated am_alc female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_health treated am_health female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
esttab , se star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("Index" "Bed" "Tired" "Alcohol" "Health") title(Regressions lifestyle index subparts (std)) drop(female_b yearofstudy_b UIB above_median_health above_median_study above_median_happiness above_median_selfcontrol am_bed am_tired am_alc am_health) replace

*** Subparts self-control (Table A7)
gen thinkbefore_m=11-acttoofast_m
gen resist_m=11-difficultresisttemptations_m
gen dontprocrastinate_m=11-procrastinate_m

sum  thinkbefore_m resist_m dontprocrastinate_m sticktoplans_m

gen std_thinkbefore_m=thinkbefore_m/2.154883
gen std_resist_m=resist_m/2.271482 
gen std_dontprocrastinate_m=dontprocrastinate_m/2.485183
gen std_sticktoplans_m=sticktoplans_m/1.926332

*above median subindex
tabstat thinkbefore_b resist_b dontprocrastinate_b sticktoplans_b, stat(p50)

gen am_think=0
replace am_think=1 if thinkbefore_b>7
replace am_think=. if thinkbefore_b==.

gen am_resist=0
replace am_resist=1 if resist_b>6
replace am_resist=. if resist_b==.

gen am_dontpro=0
replace am_dontpro=1 if dontprocrastinate_b>4
replace am_dontpro=. if dontprocrastinate_b==.

gen am_stick=0
replace am_stick=1 if sticktoplans_b>7
replace am_stick=. if sticktoplans_b==.

eststo clear
eststo: quietly reg std_selfcontrol treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_dontprocrastinate_m treated  am_dontpro female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_thinkbefore_m treated am_think female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_resist_m treated am_resist female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_sticktoplans_m treated  am_stick female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
esttab , se star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("Index" "Dont procras"  "Think before"  "Resist" "Stick plan") title(Regressions selfcontrol index subparts (std)) ///
drop(female_b yearofstudy_b UIB above_median_health above_median_study above_median_happiness above_median_selfcontrol am_think am_resist am_dontpro am_stick) replace

*** Subparts happiness (Table A8)
sum happy_general_m happy_academic_m happy_health_m happy_social_m happy_economy_m
gen std_happy_g=happy_general_m/1.944128 
gen std_happy_acd=happy_academic_m/2.138192 
gen std_happy_health=happy_health_m/2.172137 
gen std_happy_soc=happy_social_m/2.117654
gen std_happy_econ=happy_economy_m/2.434875 

eststo clear
eststo: quietly reg std_happy_g treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_happy_acd treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_happy_health treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_happy_soc treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg std_happy_econ treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
esttab , se stat(r2 N) star(* 0.1 ** 0.05 *** 0.01) drop(female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol) nonumbers mtitles("General" "Academic" "Health" "Social" "Economy") title(Regressions: Happiness (std)) replace

*** Correlation CSP vs baseline Lifestyle/Self-control (Table A9)
sum selfcontrol_b health_index
gen std_lifestyle_b=health_index/2.081653
gen std_selfcontrol_b=selfcontrol_b/ 1.48796

eststo clear
eststo: quietly reg std_CSP std_lifestyle_b female_b yearofstudy_b UIB std_RSP if treated==0,r
eststo: quietly reg std_CSP std_selfcontrol_b female_b yearofstudy_b UIB std_RSP if treated==0,r
eststo: quietly reg std_CSP_F15 std_lifestyle_b female_b yearofstudy_b UIB std_RSP_F15,r
eststo: quietly reg std_CSP_F15 std_selfcontrol_b female_b yearofstudy_b UIB std_RSP_F15,r

esttab , se star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("stdCSP" "stdCSP" "stdCSP" "stdCSP") drop(female_b yearofstudy_b UIB std_RSP std_RSP_F15)

*** Hetero academic controlling for F15 (Table A10)
eststo clear
eststo: quietly reg std_CSP treated above t_above female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_CSP std_CSP_F15 treated above t_above female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_CSP std_GA_F15 std_CSP_F15 treated above t_above female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)

eststo: quietly reg std_GA treated above t_above female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_GA std_GA_F15 treated above t_above female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_GA std_CSP_F15 std_GA_F15 treated above t_above female_b yearofstudy_b UIB std_RSP, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)

esttab , se stat(skal skalse N) star(* 0.1 ** 0.05 *** 0.01) drop(female_b yearofstudy_b UIB std_RSP) nonumbers mtitles("CSP" "CSP" "CSP" "CSP" "CSP" "GA") title(Regressions: heterogeneity (std)) replace

*** Intermediary variables Hetero (Table A11)
eststo clear
eststo: quietly reg std_lifestyle treated t_above_health above_median_study  above_median_health  above_median_happiness  above_median_selfcontrol female_b yearofstudy_b UIB, r
lincom treated+t_above_health
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_selfcontrol treated t_above_self above_median_study  above_median_health  above_median_happiness  above_median_selfcontrol female_b yearofstudy_b UIB, r
lincom treated+t_above_self
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_happy treated t_above_happy above_median_study  above_median_health  above_median_happiness  above_median_selfcontrol female_b yearofstudy_b UIB, r
lincom treated+t_above_happy
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_study treated t_above_study above_median_study  above_median_health  above_median_happiness  above_median_selfcontrol female_b yearofstudy_b UIB, r
lincom treated+t_above_study
estadd scalar skal r(estimate)
estadd scalar skalse r(se)

esttab , se stat(skal skalse N) star(* 0.1 ** 0.05 *** 0.01) drop(female_b yearofstudy_b UIB) nonumbers mtitles("Lifestyle" "Self-control" "Happiness" "Study hours") title(Regressions: heterogeneity (std)) replace

*** Gym attendance hetero (Table A12)
eststo clear
eststo: quietly reg std_gymvisits treated above_median_health t_above_health above_median_study above_median_happiness above_median_selfcontrol female_b yearofstudy_b UIB, r
lincom treated+t_above_health
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_gymvisits treated above_median_selfcontrol above_median_study above_median_health above_median_happiness t_above_self female_b yearofstudy_b UIB, r
lincom treated+t_above_self
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_gymvisits treated above_median_happiness above_median_study above_median_health  above_median_selfcontrol t_above_happy female_b yearofstudy_b UIB, r
lincom treated+t_above_happy
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_gymvisits treated above_median_study t_above_study above_median_health above_median_happiness above_median_selfcontrol female_b yearofstudy_b UIB, r
lincom treated+t_above_study
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_gymvisits treated above t_above female_b yearofstudy_b UIB, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)

esttab , se stat(skal skalse N) star(* 0.1 ** 0.05 *** 0.01) drop(female_b yearofstudy_b UIB) nonumbers mtitles("Gym" "Gym" "Gym" "Gym" "Gym" "Gym") title(Regressions: heterogeneity (std)) replace

*** Academic performance by subtreatment (Table A13)
eststo clear
eststo: quietly reg std_gymvisits card pt bonus female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
test card=pt
estadd scalar p_diff1=r(p)
test card=bonus
estadd scalar p_diff2=r(p)
test pt=bonus
estadd scalar p_diff3=r(p)
eststo: quietly reg std_CSP card pt bonus female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP,r
test card=pt
estadd scalar p_diff1=r(p)
test card=bonus
estadd scalar p_diff2=r(p)
test pt=bonus
estadd scalar p_diff3=r(p)
eststo: quietly reg std_GA card pt bonus female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP,r
test card=pt
estadd scalar p_diff1=r(p)
test card=bonus
estadd scalar p_diff2=r(p)
test pt=bonus
estadd scalar p_diff3=r(p)
esttab , se star(* 0.1 ** 0.05 *** 0.01) stat(p_diff1 p_diff2 p_diff3)  nonumbers mtitles("Gym" "CSP" "GA" ) title(Regressions: subtreat) drop(female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP) replace

***CSP Retake and Pass/Fail (Table A14)
eststo clear
eststo: quietly reg std_CSP_retake treated, r
eststo: quietly reg std_CSP_retake treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP, r
eststo: quietly reg std_CSP_pf treated, r
eststo: quietly reg std_CSP_pf treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP, r
esttab , se stat(N) star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("CSP (retake)" "CSP (retake)" "CSP (P/F)" "CSP (P/F)") title(Regressions: Robustness (std)) replace

*** Self-reported exercise (Table A15)
gen exercise_dummy=0 if timesexercised_m==0
replace exercise_dummy=1 if timesexercised_m!=0 & timesexercised_m!=. 

eststo clear
eststo: quietly reg std_exercise treated,r
eststo: quietly reg std_exercise treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
eststo: quietly reg exercise_dummy treated,r
eststo: quietly reg exercise_dummy treated female_b yearofstudy_b UIB above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
esttab , se star(* 0.1 ** 0.05 *** 0.01) nonumbers mtitles("Times Exercise (std)" "Times Exercise (std)" "Exercised (1/0)" "Exercised (1/0)") title(Regressions: Gym (std)) replace

*** Hetero self-reported exercise (Table A16)  
eststo clear
eststo: quietly reg std_exercise treated above_median_health t_above_health  above_median_selfcontrol above_median_study above_median_happiness female_b yearofstudy_b UIB, r
lincom treated+t_above_health
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_exercise treated above_median_selfcontrol t_above_self above_median_health  above_median_study above_median_happiness female_b yearofstudy_b UIB, r
lincom treated+t_above_self
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_exercise treated above_median_happiness t_above_happy above_median_health above_median_selfcontrol above_median_study female_b yearofstudy_b UIB, r
lincom treated+t_above_happy
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_exercise treated above_median_study t_above_study above_median_health above_median_selfcontrol above_median_happiness female_b yearofstudy_b UIB, r
lincom treated+t_above_study
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_exercise treated above t_above female_b yearofstudy_b UIB, r
lincom treated+t_above
estadd scalar skal r(estimate)
estadd scalar skalse r(se)

esttab , se stat(skal skalse N) star(* 0.1 ** 0.05 *** 0.01) drop(female_b yearofstudy_b UIB) nonumbers mtitles("Gym" "Gym" "Gym" "Gym" "Gym" "Gym") title(Regressions: heterogeneity (std)) replace

*** Effect by studentpopulation (Table A17)
gen t_uib=treated*UIB

eststo clear
eststo: quietly reg std_CSP treated UIB t_uib std_RSP female_b yearofstudy_b above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
lincom treated+t_uib
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
eststo: quietly reg std_GA treated UIB t_uib std_RSP female_b yearofstudy_b above_median_study above_median_health above_median_happiness above_median_selfcontrol,r
lincom treated+t_uib
estadd scalar skal r(estimate)
estadd scalar skalse r(se)
esttab , se star(* 0.1 ** 0.05 *** 0.01)  stat(skal skalse N)  nonumbers mtitles("CSP" "GA") drop(female_b yearofstudy_b above_median_study above_median_health above_median_happiness above_median_selfcontrol std_RSP) title(Regressions: Academic by student(std)) replace


*******************
*** Figures in paper/appendix
*******************

*******
*Paper
*******

***Cumulative Gym (Figure 2)
gen vps_c=visitspersemester
replace vps_c=20 if vps_c>20 & vps_c!=.

cumul vps_c if treated==0, gen(cumul_vps_c)
cumul vps_c if treated==1, gen(cumul_vps_t)

twoway (line cumul_vps_c vps_c, sort(cumul_vps_c cumul_vps_t vps_c) ///
connect(J) lpattern(l) lcolor(black)) (line cumul_vps_t vps_c, ///
sort(cumul_vps_c cumul_vps_t vps_c) connect(J)  lpattern(_) lcolor(black)), ytitle(Cumulative distribution) xtitle(Gym visits) legend(label(1 "Control") label(2 "Treatment"))

graph export Fig2.eps, replace

*** Barrier to exercise (Figure 3)

/*produce estimates for figure and save beta/se's in excelsheet "datafigure3.xlsx"
reg std_gymvisits treated if expensivebuttime_8_na==0, r 
reg std_gymvisits treated if expensivebuttime_8_na==1, r 

reg std_gymvisits treated if notime_8_na==0, r 
reg std_gymvisits treated if notime_8_na==1, r 
*/

preserve
import excel "./datafigure3.xlsx", firstrow clear

destring coeff se, replace
gen lo=coeff-se
gen hi=coeff+se

gen group2=group
replace group2=. if group2==2
replace group2=2 if group2==4
replace group2=. if group2==3
 
twoway (scatter coeff group2 if inlist(group,1,4)) ///
	(rcap lo hi group2 if inlist(group,1,4)), ///
	ylabel(0(0.2)1.2, nogrid) xlabel(0.5 " " 1 "Expensive but time" 2 "No time" 2.5 " ", noticks)  ///
	xtitle("Stated barrier to exercise") ytitle("Treatment effect") subtitle(Gym visits (std)) ///
	yline(0.629) legend(label(1 "Coefficient") label(2 "Robust S.E."))

graph export Fig3.eps, replace

***Treatment effect Academic (Figure 4)	

/*produce estimates for figure and save beta/se's in excel sheet "datafigure4.xlsx"
reg std_CSP treated,r
reg std_GA treated,r
*/

preserve
import excel "datafigure4.xlsx", firstrow clear

destring coeff se, replace
gen lo=coeff-se
gen hi=coeff+se

twoway (scatter coeff GA) ///
	(rcap lo hi GA), ///
	ylabel(-0.3(0.1)0.3, nogrid) xlabel(-0.5 " " 0 "CSP (std)" 1 "GA (std)" 1.5 " ", noticks)  ///
	xtitle("") ytitle("Treatment effect") ///
	yline(0) legend(label(1 "Coefficient") label(2 "Robust S.E."))

graph export Fig4.eps, replace

***Hetro treatment effect Academic (Figure 5)

/*produce estimates for figure and save beta/se's in excel sheet "datafigure5.xlsx"
reg std_CSP treated if above==0,r
reg std_CSP treated if above==1,r

reg std_GA treated if above==0,r
reg std_GA treated if above==1,r
*/
* net install grc1leg2.pkg, from (http://digital.cgdev.org/doc/stata/MO/Misc/)
preserve
import excel "datafigure5.xlsx", firstrow clear

destring coeff se, replace
gen lo=coeff-se
gen hi=coeff+se

twoway (scatter coeff above if outcome=="CSP") ///
	(rcap lo hi above if outcome=="CSP"), ///
	ylabel(-1(0.2)1, nogrid) xlabel(-0.5 " " 0 "Above median: No " 1 "Above median: Yes" 1.5 " ", noticks)  ///
	xtitle("") ytitle("Treatment effect") subtitle(Panel A: Completed Study Points (std)) ///
	yline(0) legend(label(1 "Coefficient") label(2 "Robust S.E."))saving(csp.gph)

twoway (scatter coeff above if outcome=="GA") ///
	(rcap lo hi above if outcome=="GA"), ///
	ylabel(-1(0.2)1, nogrid) xlabel(-0.5 " " 0 "Above median: No " 1 "Above median: Yes" 1.5 " ", noticks)  ///
	xtitle("") ytitle("Treatment effect") subtitle(Panel B: Grade Average (std)) ///
	yline(0) legend(label(1 "Coefficient") label(2 "Robust S.E.")) saving(ga.gph)
 	
grc1leg2 csp.gph ga.gph

graph export Fig5.eps, replace

******************
*Appendix Figures
******************

*** Barrier to exercise full scale (Figure A1)
preserve
collapse (mean) m = std_gymvisits (sem) se=std_gymvisits, by(treated SIBcardexpensive_b)
gen hi=m+se
gen lo=m-se
twoway (scatter m SIBcardexpensive_b if treated==0) (scatter m SIBcardexpensive_b if treated==1) ///
(rcap hi lo SIBcardexpensive_b if treated==0) (rcap hi lo SIBcardexpensive_b if treated==1), ///
legend(order(1 2)) legend(label(1 "Control")) legend(label(2 "Treated")) ///
xtitle("SIB card expensive (1-10)") xlabel(1(1)10,) ytitle(Gym visits (std)) subtitle(Panel A: Cost Barrier) saving(cost, replace) 

preserve
collapse (mean) m = std_gymvisits (sem) se=std_gymvisits, by(treated notimetoexercise_b)
gen hi=m+se
gen lo=m-se
twoway (scatter m notimetoexercise_b if treated==0) (scatter m notimetoexercise_b if treated==1) ///
(rcap hi lo notimetoexercise_b if treated==0) (rcap hi lo notimetoexercise_b if treated==1), ///
legend(order(1 2)) legend(label(1 "Control")) legend(label(2 "Treated")) ///
xtitle("No time to exercise (1-10)") xlabel(1(1)10,) ytitle(Gym visits (std))  subtitle(Panel B: Time Barrier) saving(time, replace) 

graph combine cost.gph  time.gph, col(1)

graph export FigA1.eps, replace

***Cumulative Academic (Figure A2)
gen csp_c=CSP
replace csp_c=40 if csp_c>40 & csp_c!=.

cumul csp_c if treated==0, gen(cumul_csp_c)
cumul csp_c if treated==1, gen(cumul_csp_t)

twoway (line cumul_csp_c csp_c, sort(cumul_csp_c cumul_csp_t csp_c) ///
lpattern(l) connect(J) lcolor(black)) ///
(line cumul_csp_t csp_c, sort(cumul_csp_c cumul_csp_t csp_c) ///
connect(J)  lpattern(_) lcolor(black)), ///
ytitle(Cumulative distribution) xtitle(CSP) legend(label(1 "Control") label(2 "Treatment")) subtitle(Panel A) saving(CSP, replace) 

cumul GA if treated==0, gen(cumul_ga_c)
cumul GA if treated==1, gen(cumul_ga_t)

twoway (line cumul_ga_c GA, sort(cumul_ga_c cumul_ga_t GA) ///
lpattern(l) connect(J) lcolor(black)) ///
(line cumul_ga_t GA, ///
sort(cumul_ga_c cumul_ga_t GA) connect(J)  lpattern(_) lcolor(black)), ///
ytitle(Cumulative distribution) xtitle(GA) legend(label(1 "Control") label(2 "Treatment")) subtitle(Panel B) saving(GA, replace) 

grc1leg2 CSP.gph  GA.gph

graph export FigA2.eps, replace


*** Gym by week (Figure A3)
preserve
use gymvisits_perweek.dta,clear

twoway (connected mean_gymvisits week_nr if treated==0) ///
	(connected mean_gymvisits week_nr if treated==1), ///
	ylabel(0(0.2)0.9,) xlabel(1(5)25,) text(0.87 3.5 "Pre-intervention") text(0.87 12 "Teaching period") ///
	text(0.87 22 "Exam period") xtitle("Week of Year") ///
	xline(7, lcolor(black)) xline(18, lcolor(black)) ///
	ytitle(Visits per Week) legend(label(1 "Control") label(2 "Treatment")) subtitle(Gym Visits per Week)

graph export FigureA3.eps,replace
